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ABSTRACT 


One of NASA' 3 primary mission objectives for its first* flights of the 
Space Shuttle System is to revisit the Sky lab vehicle. The purpose of this 
revisit will be to add a kick stage to the vehicle to raise its orbit and to 
do a limited number of experiments. Without the addition of the kick motor 
Skylab's orbit will decay by mid-1980. Unfortunately, it is highly likely 
that the Sky lab will be in a tumbling mode, because it is not controlled and 
it is under the influence of environmental effects. The work reported in the 
proposed paper deals vith the prediction of this uncontrolled motion. Com- 
puter simulations were employed to generate results based on a fairly com- 
plicated analytical model. Environmental effects considered Include atmospheric 
drag and gravity gradients. Mathematical models were based on Euler's moment 
equations and transformations from the body frame to inertial frame. A range 
of altitudes was considered, from 278 kilometers to 417 kilometers. Sky lab 
is predicted to be in this orbital range during the 1979-1980 period. Results 
indicated that there will be a cyclic attitude motion due to the combined ef- 
fect of gravity gradient and atmosphere drag. This motion can be characterized 
as being a slow tumble. A three-degree-of-freedom simulator is used to show 
what the space shuttle crew will see as they approach the Skylab. 
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Dynamic pressure - 1/2 p v 
Satellite distance from earth's center 
Satellite orbital velocity 

Satellite velocity relative to the incident stream 
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Satellite angle of attack 

Earth's gravitational constant 
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Atmospheric density 

Density cosine curve fit constants 
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Satellite angular velocity 
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CHAPTER 1 


Introduction 

1.1 Origin and Importance of Investigation 

A revisit to the Skylab spacecraft is one of the early priority missions 
NASA envisions for the Space Shuttle. Currently Skylab is in a decaying orbit 
around the earth such that by late 1979-1980 it will enter the atmosphere and 
burn up. The purpose of the Space Shuttle visit would be to add a kick stage 
to Skylab to raise its orbit so it could be used for further experimentation. 
However, Skylab is now passive and is under the influence of environmental 
perturbations which cause it to tumble. Skylab must be motionless in order 
that work may be done on the vehicle. Therefore, an idea of what Skylab will 
look like as the Space Shuttle crew approaches is necessary so that proper 
and adequate equipment to slow or stop the tumble can be placed on board the 
Shuttle. 

1.2 Statement of the Problem and Scope of Investigation 

The purpose of this investigation, then, is to predict Sky lab’s attitude 
motion in 1979-1980. Various perturbation forces acting on the passive ve- 
hicle are considered and numerically analyzed. These expressions are incor- 
porated into Euler's moment equations and then solved via a digital computer. 
The moment equations are solved for four altitudes ranging from 278 kilometers 
to 417 kilometers since orbital decay prediction for Skylab shows this range 
of altitude for Skylab* s orbit at the time of the Shuttle visit. The resulting 
solutions show cyclic tumbling. A three-degree-of-freedom simulator is used 
to visualize the predicted motion of Skylab. 
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This investigation lays groundwork for NASA to decide If a Skylab revisit 
is feasible. It also provides information that can be of great Importance 
for predicting other satellite tumbling modes so that they, too, may be re- 
visited and reactivated. 
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CHAPTER 2 

Background Information 

2.1 Coordinate Systems 

In order to begin a study of the motion of a satellite, a coordinate 
system to describe the position and attitude of the satellite must be established. 
Two such systems are here defined. 

2.1.1 Inertial Reference Coordinate System 

Figure 1 depicts the standard geocentric inertial coordinate system. It 
has its origin at the earth’s center with the Xj - Y^ plane coinciding with 
the earth's equatorial plane and Xj pointing toward the Vernal Equinox. The 
inertial reference coordinate is denoted in matrix form 

2.1.2 Satellite Body Coordinate System 
Body coordinates, 


are related to inertial coordinates through the Euler angles i{>, 0, as shown 
in figure 2. The transformation matrix for inertial to body coordinated is 

(cos <J> cos if) - sin $ cos 0 sin tJO 
(-sin 4> cos ip - cos cos 0 sin i|0 
(sin 0 sin ij>) 


derived in reference (1) as: 

X 
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( 2 . 1 . 2 - 1 ) 

Of special Interest to this investigation Is the transformation that relates 
the Euler rates to a body angular velocity vector, u>, also shown In Figure 2. 
This transformation matrix is found in reference (1) as: 
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2.2 Equations of Motion 

In reference (1), Euler's equations of motion for principal axes are given 
as: 

M - A u + 03 0) (C-B) (2. 2-la) 

x x y z 7 

M - B ai +0) 0) (A-C) (2.2-lb) 

y y x z 

M - C to + 0 ) 0 ) (B-A) (2.2-lc) 

z z x y 

Where A, B, and C are the principal moments of inertia coinciding with the 
principal body axes x, y, and z, u is the angular velocity, and w is the 
angular acceleration. It is important to note that equations (2.2-1) are 
expressed in body coordinates. 
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CHAPTER 3 

Discussion, Analysis, and Reduction of 
External Moments 

Reduction of environmental perturbation moment expressions to a manageable 
form requires theoretical analysis to determine which outside forces are small 
and can be neglected without loss of accuracy. The following deals with dis- 
cussion and reduction of the four principal perturbing forces — gravity gradient, 

V 

atmospheric drag, solar radiation, and magnetic torques — so that these forces 
can be incorporated Into a digital program to determine Euler angles and rates. 

3.1 Gravity Gradient 

The gravitational attraction between the enormous mass of the earth and 
a satellite in a fairly close orbit about the earth produces a torque which 
is a major perturbation of the attitude of the satellite. The gradient of the 
gravity force produces this perturbation. The gravity gradient is a function 
of satellite altitude, satellite orientation, earth mass distribution, and 
satellite mass distribution. 

By assuming a spherical earth of relatively uniform mass distribution, 

and by assuming that the satellite is a point mass as compared to the mass of 

the earth, expressions for the moments produced from the gravity gradients 

are illustrated by Tate in reference (2) as follows: 

# 

M - sin 20 cos 2 9 (C-B) (3.1-la) 

831 2R* 
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M - sin 26 cos <J> (C-A) (3.1-lb) 

* 2R 2 


M - -^~r sin 26 sin $ (A-B) (3.1-lc) 

gZ 2R 2 


One important note is that equations (3.1-1) are derived and written in body 
coordinates, and R is the distance from the center of the earth (the origin 
of the geocentric coordinate system) to the satellite center of mass. Values 
of A, B, and C, the satellite principal moments of inertia, are given in 
appendix A. 


3.2 Aerodynamic Drag 

Around the year 1980, Skylab will be orbiting at an altitude of less than 
350 kilometers. Aerodynamic drag, a function of atmospheric density, satellite 
velocity, satellite shape, and satellite angle of attack, becomes a major ex- 
ternal perturbation force at altitudes of 800 kilometers or less. Therefore, 
in this investigation, aerodynamic drag moments must be considered. Following 
the format in reference (2) , external moments due to aerodynamic drag are found 
to be: 
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where c is the roll moment coefficient, c is the pitch moment coefficient, 
x 7 

c z is the yaw moment coefficient, A^^ is the Sky lab reference area, is 

1 2 

the Skylab reference diameter, q is the dynamic pressure, -? p V q , and V q is 
the orbital velocity. 

Aerodynamic drag depends upon atmospheric density, p, and the pi ch, roll, 
and yaw moment coefficients. However, density is a function not only of al- 
titude but also of other variations that change with time. The moment coef- 
ficients vary with time due to Skylab's changing attitude. In order to 
accomodate these time dependent variations, density and moment coefficients 
must be handled separately and then incorporated into the atmospheric drag 
moment expressions. 

Interpolation of information from the chart of Figure 3 indicates that 
Skylab will be within 278 to 417 kilometers above the earth during the time 
in question (1979-1980). Therefore, four representative altitudes of 278, 

324, 370, and 417 kilometers were chosen. Evaluations of all perturbations were 
made at these four altitudes in order to produce a more general overview of 
the probable motion of Skylab. 

3.2.1 Dens ity Evaluation 

Since the possible orbital altitudes of Skylab are well over 90 kilometers, 
the Jacchia (1970) model for atmospheric density at satellite altitudes was 
chosen to evaluate the density variations of the orbits. The Jacchia model 
in its completeness accounts for temperature and density variations due to 
seasonal and latitudinal variations, to both solar and geomagnetic activity, 
and diurnal and semi-annual variations. A copy of the basic Jacchia program 
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is round In reference 4. The Jacchia model program as supplied by NASA had 
co be slightly reworked and adapted to fit the constraints of the IBM 370/168 
computer used at The Pennsylvania State University. 

The Jacchia model will calculate density values •w" minute, every day. 

In order to reduce computer time, noon March, 1980 was chosen as an average 
reference time to represent early 1980. The solar indices for this date were 
obtained from NASA. The density distribution about the earth was then cal- 
culated using the Jacchia program. This distribution was modelled to represent 
a cosine curve illustrating the higher density values on the sunside of the 
earth. These cosine curve fits generally vrere accurate to within five percent 
of the actual calculated values. The density distribution about the earth 
could then be obtained from the following equation 

p - p_ + p. cos (a) t) (3, 2. 1-1) 

BA O 

The values of p^, and confer each reference altitude are given in appendix 
A. The variable t in equation (3. 2. 1-1) is time, where the time relates to 
the longitudinal position in the orbit. 


3.2.2 Moment Coefficient Evaluation 

A Fourier series curve fit formula, obtained from NASA by Tate (2), cal- 
culates the aerodynamic drag moment coefficients as functions of the satellite 

angle of attack, ct , and the roll angle, <p , as shown in Figure 4. 

A A 

The drag coefficient equation is: 


A m 

C(ct a , ^a ) " T + E * A i COS 1 a a + B i Sin 1 °a^ 


(3.2. 2-1) 
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A 1 <V ' ^ + ? KlJ cos 1 ♦. + b aij 9ln 3 V < 3 - 2 - 2 ' 2 > 

J-l 
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+ ^! C0S 3 *« + b biJ 8i ” 3 '‘a 1 


(3.2. 2-3) 


Fourier series analysis leads to the finding that coefficients evaluated 
at 1 * 1,3 and j ■ 2,4,6 and coefficients evaluated at 1 * 1,3, and j ■ 1,3,5, 
yield a highly accurate representation of the coefficient curve. 

The a^ A^, b^A^, a ij B i and b ij B i coe ^^ icient values are given in Tate 
(2). These coefficient values are based on Sky lab with the auxiliary thermal 
shield, ATM solar arrays, and orbital workshop solar panel No. 1 deployed as 
shown in Figure 4. 

The angles $ and ct a are ev^ luated xn body coordinates such that 

0° < (J> < 360° (3. 2. 2-4) 
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(3.2. 2-5) 


Since the orbit of Sky lab is assumed circular, V 


V where V is given 
o o 


in body coordinates by 
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V ■ V cos 8 cos 'P (3.2.2-6a) 

X o 

V “V sin <J> sin 0 cos $ - cos $ sin V (3.2.2-6b) 

y o 

V z ■ V q cos <J> sin 0 cos + sin <p sin (3.2.2-6c) 


The equations for calculating the aerodynamic moment coefficients, c , c , 

x y 

and c , were written into a separate subroutine to be used in the subroutine 
z 

to calculate total aerodynamic drag moments 


3.2.3 Complete Subroutine for Aerodynamic Drag 

The first part of the subroutine to calculate aerodynamic drag at ^he 

reference altitudes incorporated the cosine curve fit equations as a function 

of height and time. The second part consisted of the subroutine to calculate 

the moment coefficients as a function of velocity, angle of attack, and roll 

angle, as well as the reference area, A rgf , and the reference dicmeter, & ref » 

as given in Appendix A. By incorporating equations (3.2-1) and the expression 

1 2 

for aerodynamic pressure, q • -j p V q , the expressions for aerodynamic moments 
are obtained as 
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in body coordinates. 



3.3 Solar Wind Radiation 

The sun produces radiation which will strike a body in space and impart 
momentum to that body. The amount of momentum transferred to the body in the 
form of torques depends on the distance the body is from the sun, the body 
surface geometery, and the body surface reflectivity. Examination of the moments 
due to solar wind and solar radiation indicate that their effect on Sky lab will 
remain negligible. Since solar effects are relatively small at low altitudes 
and since the solar moment expression is difficult to handle numerically, the 
effect of solar moments can be neglected while still achieving a good deal of 
accuracy in the modelling of Sky lab's attitude motion. 

3.4 Magnetic Torques 

A satellite's magnetic components interact with the earth's magnetic 
field causing perturbation torques on the satellite. However, since Sky lab 
is passibe, the magnetic interaction between it and the earth are assumed 
negligible, and are not included in this investigation. 

3.5 Final Moment Equations 

The final moment equations that were numerically integrated are found 
by combining equations (3.1-1) and 0.2. 3—1) with equation (2.2-1) such that 

M +M * A u> + u to (C-B) (3.5-la) 

gx ax x y z 

M + M * A <u + w a) (A-C) (3.5-lb) 

gy ay y x z 

M + M = A <o + <o co (B-A) 

gz az z x y 


(3.5-lc) 
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CHAPTER 4 
Results 

Euler's equations of motion, equations (3.5-1), were solved numerically 
using an IBM 370/168 digital computer. The external moments were calculated 
by the methods described in Chapter 3. The integration routine consisted of 
a second order Runge-Kutta method used for the first four iterations to start 
the solution. The routine then switched to a Hamming predictor-corrector 
method for the remainder of the program. A time increment of five seconds 
was used. 

The input for the program was the initial Euler angles, the orbital 
altitude, the density model and constants, and the Fourier aerodynamic moment 
model and constants. The input can be found in appendix A. The output of 
the integration routine was the body rates and their time derivatives as 
functions of time. Equation (2.1. 2-2) was then used to transform the body 
rates to the final output, the Euler rates. Using a more standard notation 
the body, x, y, and z rates may be respectively denoted as the yaw, roll, and 
pitch rates with respect to the orbital plane as shown in Figure 5. 

Figures 6 through 9 are plots of the precession, nutation and spin rates 
as functions of time for each of the four altitudes considered. After one 
orbital period, all precession and spinning motion show cyclic tumbling 
characteristic of gravity gradient perturbation. The notation motion lacks 
this cyclic behavior because there is no gradient of gravitational force in 
this mode. This can be seen in the relatively slow increase and small nutation 
rate which implies only aerodynamic torques affect this mode. 
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The singular points In the precession and spin rates are caused by 6 the 
nutation angle, approaching zero. The transformation used, equation (2. 1.2-2), 
has a factor of — - — r which causes singularity for 0 approaching zero. In 

SIB 0 

the case of the 324 kilometer plots, 6 actually passes through zero, delaying 
the cyclic tumbling in the precession and spin rates. 

Examination of the plots for 278 kilometers. Figure 6, shows that the 
precession and spin rates reach a steady state tumbling motion after about one 
orbital period, approximately 92 minutes. Ignoring the singular points, the 
rates and periods of this steady state motion are: 

ip. Precession rate: Period » 40 minutes 

Rates * -.8 to .8 degrees/ second 
0, Nutation rate: Motion is not cyclic 

Rates ■ -.4 to .4 degrees/second 
<p , Spin rate: Period - 40 minutes 

Rates * -.5 to .5 degrees /second 

A three degree of freedom simulator was rsed to better visualize the at- 
titude motion. The Euler rates obtained from the digital computer were trans- 
formed to DC voltage using a digital to analog converter. These voltages then 
fed the motors which ran the gimbals of the simulator. The configuration is 
shown in Figure 10. 
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APPENDIX A 


Satellite Parameters 


Principal Moments of Inertia 

5 2 

A * 7.93321 x 10 kilograms -meter 

6 2 

B “ 3.767828 x 10 kilograms - meter 

6 2 

C * 3.694680 x 10 kilograms - meter 

Aerodynamic Reference Area and Diameter 

2 

A 79.46 meter 
ref 

2 

D r = 10.058 meter 
ref 

Density Reference Pete and Time 
March 1, 1-980 
12:00 noon 

Julian Day = 2444301 
Density Constants 

417 kilometers u) =* .06147591 seconds ^ 
o 

-12 3 

p « 2.785 x 10 kilogram/meter 
8 

‘-12 3 

p ■ 1.295 x 10 kilogram/ meter 

A 

370 kilometers w =* .0637886833 seconds ^ 
o 

12 3 

p_ ■ 6.3 x 10 kilogram/meter 

D 

12 3 

p^ ■ 2.55 x 10 kilogram/ meter 



APPENDIX A (continued) 


324 kilometers w 

o 

P B 

P A 

287 kilometers ^ 

o 

P B 


.0646639997 seconds 
1.4975 x 10 kilogram/meter 3 
5.135 x 10*2 kiiogram/meter 3 
.0654498469 seconds * 

3.9295 x 10 ** kilogram/meter 3 
1.0785 x 10”** kilogram/meter 3 


Initial Euler Angles 
^ ® 1.07 degrees 
0 * -79.96 degrees 
0 * 12.85 degrees 




to Vernal Equatorial Plane 

Equinox 


Figure 1. Geocentric Inertial Coordinate System 
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Figure 2. Euler angles, rates, and body rates 
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Figure 5. Pitch, 
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Figure 6. Precession, Nutation, and Spin Rates for 278 Kilometers 
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Figure 10. Simulator with Euler rates 




